function [ymin thetamin init lower upper] = launch_optimization(optimizer)

% Note that ordering here must match likelihoods.cpp

eps = 1e-5; % TODO: figure out where numeric instability in C code is
    
% Initial parameters
init = [ ...
    0.8 ...  % wall penalty
    0.6 ...  % occlusion penalty
    100 ...  % "ideal" features
    100 ...
    100 ...
    100 ...
    0.1 ...    % "fall-offs" (w.r.t normalized y-coordinates in [0,1])
    0.1 ...
    0.1 ...
    0.1 ...
    10 ...   % noise bandwidth
    10 ...
    10 ...
    10 ...
    0.5 ...  % foreground probabilities
    0.5 ...
    0.5 ...
    0.5 ...
    ];

% Lower bounds
lower = [ ...
    eps ...  % wall penalty
    eps ...  % occlusion penalty
    0 ...    % "ideal" features
    0 ...
    0 ...
    0 ...
    1e-3 ...  % "fall-offs" (w.r.t normalized y-coordinates in [0,1])
    1e-3 ...
    1e-3 ...
    1e-3 ...
    1 ...    % noise bandwidth
    1 ...
    1 ...
    1 ...
    eps ...  % foreground probabilities
    eps ...
    eps ...
    eps ...
    ];

% Upper bounds
upper = [ ...
    1-eps ...  % wall penalty
    1-eps ...  % occlusion penalty
    1000 ...   % "ideal" features
    1000 ...
    1000 ...
    1000 ...
    2 ...    % "fall-offs" (w.r.t normalized y-coordinates in [0,1])
    2 ...
    2 ...
    2 ...
    500 ...    % noise bandwidth
    500 ...
    500 ...
    500 ...
    1-eps ...  % foreground probabilities
    1-eps ...
    1-eps ...
    1-eps ...
    ];

% Transform between log-likelihood and objective function for optimization
% Log-likelihood is always <= 0
%ty = @(y)(-y);
inv_ty = @(y)(-y);


% Report objective function value at initial point
yinit = f(init);
disp(['Objective function at init: ' num2str(yinit)]);

if (strcmp(optimizer, 'direct'))
    disp('Optimizing with Direct...');
    Problem.f = @f;
    opts.maxevals = 1000;
    opts.showits = true;
    [ymin, thetamin] = Direct(Problem, [lower' upper'], opts);

elseif (strcmp(optimizer, 'gpgo'))
    disp('Optimizing with gpgo...');
    opt.save_str = 'opt_alex_flint';
    opt.function_evaluations = 500;
    opt.exp_loss_evals = 100;
    opt.num_hypersamples = 1;
    opt.derivative_observations = false;
    [ymin thetamin] = gpgo(@f, init, lower, upper, opt);
else
    error(['Unknown optimizer: ' optimizer]);
end

% Flip back
ymin = inv_ty(ymin);

function [neglogL, grad_neglogL] = f(theta)
% Flip problem upside down (for minimization)
neglogL = -compute_loglik(theta);
if nargout>1
[dummy, grad_neglogL] = compute_loglik(theta);
grad_neglogL = -grad_neglogL;
end